Male courtship behaviour - aDN > TNT

Setup

library("tidyverse")
library("ggpubr")
library("zoo")
setwd("/mnt/LocalData/behaviour/aDN/aDN_behaviour")

Colours …

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}

gg_color_hue(2)

[1] "#F8766D" "#00BFC4"

genotypes <- read_tsv("../2019_03_06_Courtship/genotype.tsv",col_names = TRUE)
genotypes

# indices_list <- list.files("raw data/") %>% str_subset("_Indices.csv")
indices_list <- list.files("../2019_03_06_Courtship/",recursive = TRUE) %>% str_subset("_Indices.csv") %>% str_subset("Male")
all_indices <- tibble()
for (indices_file in indices_list) {
  video_name <- indices_file %>% str_remove("/.*")
  temp <- read_csv(paste0("../2019_03_06_Courtship/",indices_file))
  temp <- temp %>%
    mutate(video = video_name)
  all_indices <- bind_rows(all_indices,temp)
}
all_indices

all_male_indices <- left_join(x = genotypes,y = all_indices,by=c("video"="video","fly_id"="FlyId"))

ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndexWithFacing)) + 
  geom_boxplot()

p1 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndex)) + 
  geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndexWithFacing)) + 
  geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=TotalCCI)) + 
  geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=TotalCCIwFacing)) + 
  geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4),
          # labels = c("ApproachingIndex",
          #            "ContactIndex",
          #            "EncirclingIndex",
          #            "FacingIndex",
          #            "TurningIndex",
          #            "WingIndex"),
          # hjust = 1,
          ncol = 4,
          nrow = 1)

p1 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingIndex)) + 
  geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=ContactIndex)) + 
  geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingIndex)) + 
  geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=FacingIndex)) + 
  geom_boxplot()
p5 <- ggplot(all_male_indices,aes(x=genotype,y=TurningIndex)) + 
  geom_boxplot()
p6 <- ggplot(all_male_indices,aes(x=genotype,y=WingIndex)) + 
  geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4,p5,p6),
          # labels = c("ApproachingIndex",
          #            "ContactIndex",
          #            "EncirclingIndex",
          #            "FacingIndex",
          #            "TurningIndex",
          #            "WingIndex"),
          # hjust = 1,
          ncol = 3,
          nrow = 2)

p1 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipInitiation)) + 
  geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipTermination)) + 
  geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipDuration)) + 
  geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3),
          ncol = 3,
          nrow = 1)

df <- all_male_indices %>%
  select(genotype,CourtshipTermination) %>% 
  group_by(genotype) %>% 
  mutate(len=length(CourtshipTermination))
#ggplot(df,aes(x=CourtshipTermination,color=genotype)) + geom_step(aes(len=len,y=..y.. * len),stat="ecdf") 
ggplot(df,aes(x=CourtshipTermination,color=genotype)) + geom_step(aes(y=..y..),stat="ecdf")

df <- all_male_indices %>%
  select(genotype,CourtshipDuration) %>% 
  group_by(genotype) %>% 
  mutate(len=length(CourtshipDuration))
#ggplot(df,aes(x=CourtshipDuration,color=genotype)) + geom_step(aes(len=len,y=..y.. * len),stat="ecdf") 
ggplot(df,aes(x=CourtshipDuration,color=genotype)) + geom_step(aes(y=..y..),stat="ecdf")

p1 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingBoutLength)) + 
  geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=ContactBoutLength)) + 
  geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingBoutLength)) + 
  geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=FacingBoutLength)) + 
  geom_boxplot()
p5 <- ggplot(all_male_indices,aes(x=genotype,y=TurningBoutLength)) + 
  geom_boxplot()
p6 <- ggplot(all_male_indices,aes(x=genotype,y=WingBoutLength)) + 
  geom_boxplot()
p7 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingBoutInterval)) + 
  geom_boxplot()
p8 <- ggplot(all_male_indices,aes(x=genotype,y=ContactBoutInterval)) + 
  geom_boxplot()
p9 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingBoutInterval)) + 
  geom_boxplot()
p10 <- ggplot(all_male_indices,aes(x=genotype,y=FacingBoutInterval)) + 
  geom_boxplot()
p11 <- ggplot(all_male_indices,aes(x=genotype,y=TurningBoutInterval)) + 
  geom_boxplot()
p12 <- ggplot(all_male_indices,aes(x=genotype,y=WingBoutInterval)) + 
  geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4,p5,p6,
                          p7,p8,p9,p10,p11,p12),
          # labels = c("ApproachingIndex",
          #            "ContactIndex",
          #            "EncirclingIndex",
          #            "FacingIndex",
          #            "TurningIndex",
          #            "WingIndex"),
          # hjust = 1,
          ncol = 6,
          nrow = 2)

#all_rawdata <- full_join(x = all_rawdata, y = genotypes, by = c("FileName"="video","Id"="fly_id"))
all_rawdata <- all_rawdata %>%
  full_join(x = all_rawdata, y = genotypes, by = c("FileName"="video","Id"="fly_id")) %>% 
  replace_na(list(genotype = "CS female"))
# sum(is.na(all_rawdata$genotype))
# summary(all_rawdata$genotype)
# unique(all_rawdata$genotype)
all_rawdata

# # Extract density data to average
# p <- ggplot_build(test_plot)
# ggplot(as.data.frame(p$data[[1]]), aes(x,y)) + geom_line()
# ggplot(as.data.frame(p$data[[2]]), aes(x,y)) + geom_line()

D genotype distance to other

all_rawdata <- all_rawdata %>%  
  unite("unique_fly",FileName,Id, remove = FALSE) %>% 
  group_by(unique_fly) %>% 
  mutate(
    Multitasking = (Approaching + Encircling + Contact + Turning + WingGesture),
    MultitaskingWithFacing = (Approaching + Encircling + Facing + Contact + Turning + WingGesture),
    Courtship = ifelse(Multitasking>=1, 1, 0),
    CourtshipWithFacing = ifelse(MultitaskingWithFacing>=1, 1, 0),
    MultitaskingWithCopulation = (Approaching + Encircling + Contact + Turning + WingGesture + Copulation),
    MultitaskingWithCopulationWithFacing = (Approaching + Encircling + Facing + Contact + Turning + WingGesture + Copulation),
    CourtshipAndCopulation = ifelse(MultitaskingWithCopulation>=1, 1, 0),
    CourtshipAndCopulationWthFacing = ifelse(MultitaskingWithCopulationWithFacing>=1, 1, 0),
    SmoothedCourtship = ifelse((rollmean(Courtship, 150, fill = c(0,0,0), align = c("left")))>0.5, 1, 0),
    SmoothedCopulation = ifelse((rollmean(Copulation, 1250, fill = c(0,0,0), align = c("center")))>0.5, 1, 0),
    SmoothedDistToOther = ifelse((rollmean(ifelse(dist_to_other__mm > 2, 1, 0), 250, fill = c(1,1,NA), align = c("center")))>0.5, 1, 0)
)

D genotype distance to other

mean distance to other

median distance to other

mean facing angle

median facing angle

test_stat_tibble0 <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  group_by(genotype) %>% 
  filter(dist_to_other__mm > 2) %>% 
  #filter(max_wing_ang__rad > (35*pi/180)) %>%
  unite("individual", FileName:Id, remove = FALSE) %>% 
  group_by(individual) %>% 
  summarise(median_facing_angle = median(facing_angle__rad),
            genotype = unique(genotype))

aov0 <- aov(median_facing_angle~genotype,data = test_stat_tibble0)
summary(aov0)

             Df Sum Sq Mean Sq F value Pr(>F)  
genotype      3  1.036  0.3453   2.736 0.0468 *
Residuals   114 14.388  0.1262                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(aov0)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = median_facing_angle ~ genotype, data = test_stat_tibble0)

$genotype
           diff           lwr       upr     p adj
B-A -0.07051500 -0.3361762470 0.1951462 0.8999385
C-A  0.04970748 -0.1803176760 0.2797326 0.9426712
D-A  0.19276677 -0.0389146073 0.4244482 0.1380793
C-B  0.12022248 -0.1408020849 0.3812470 0.6275505
D-B  0.26328177  0.0007965104 0.5257670 0.0490063
D-C  0.14305929 -0.0832904123 0.3694090 0.3562046

median facing angle by sector

mean facing angle by sector

mean facing angle by wee sectors

counts in wee sectors

Total time in sector

Proportion time in sector

test_stat_tibble0 <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  filter(dist_to_other__mm > 2) %>% 
  filter(max_wing_ang__rad > (35*pi/180)) %>%
  group_by(unique_fly) %>% 
  summarise(counts = 100*sum(facing_angle__rad > pi/8)/length(Frame),
            genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)

$genotype
          diff        lwr       upr     p adj
B-A  6.6367127 -2.5650354 15.838461 0.2421768
C-A  7.3146600 -0.6527555 15.282076 0.0840396
D-A 16.6152204  8.5904380 24.640003 0.0000022
C-B  0.6779474 -8.3631993  9.719094 0.9973376
D-B  9.9785077  0.8867666 19.070249 0.0254776
D-C  9.3005603  1.4604520 17.140669 0.0131217

pairwise.t.test(test_stat_tibble0$counts, test_stat_tibble0$genotype)

    Pairwise comparisons using t tests with pooled SD 

data:  test_stat_tibble0$counts and test_stat_tibble0$genotype 

  A       B     C    
B 0.125   -     -    
C 0.055   0.845 -    
D 2.2e-06 0.020 0.012

P value adjustment method: holm 

Total time in sector

Proportion time in sector and within 10mm

test_stat_tibble0 <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  filter(dist_to_other__mm > 2) %>% 
  filter(dist_to_other__mm < 10) %>% 
  filter(max_wing_ang__rad > (35*pi/180)) %>%
  group_by(unique_fly) %>% 
  summarise(counts = 100*sum(facing_angle__rad > pi/8)/length(Frame),
            genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)

$genotype
          diff        lwr       upr     p adj
B-A  4.2169183 -3.2364908 11.670328 0.4558215
C-A  4.8104624 -1.6431383 11.264063 0.2159140
D-A 11.2681432  4.7680753 17.768211 0.0000889
C-B  0.5935441 -6.7297782  7.916866 0.9966463
D-B  7.0512248 -0.3130788 14.415528 0.0658244
D-C  6.4576807  0.1071987 12.808163 0.0446836

pairwise.t.test(test_stat_tibble0$counts, test_stat_tibble0$genotype)

    Pairwise comparisons using t tests with pooled SD 

data:  test_stat_tibble0$counts and test_stat_tibble0$genotype 

  A       B     C    
B 0.286   -     -    
C 0.163   0.833 -    
D 9.1e-05 0.056 0.046

P value adjustment method: holm 

Total count in inner ring

Proportion count in inner ring

test_stat_tibble0 <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  filter(facing_angle__rad <= pi/8) %>% 
  filter(max_wing_ang__rad > (25*pi/180)) %>%
  group_by(unique_fly) %>% 
  summarise(counts = 100*sum(dist_to_other__mm > 2 & dist_to_other__mm <= 3)/length(Frame),
            genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)

$genotype
          diff        lwr      upr     p adj
B-A  2.7255382  -8.488600 13.93968 0.9209524
C-A  2.5367905  -7.173071 12.24665 0.9040570
D-A 21.8215793  12.041805 31.60135 0.0000003
C-B -0.1887477 -11.207161 10.82967 0.9999678
D-B 19.0960411   8.015968 30.17611 0.0000986
D-C 19.2847888   9.730076 28.83950 0.0000040

Left vs Right wing

test_stat_tibble <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  #filter(genotype != "A") %>% 
  group_by(genotype) %>% 
  group_by(unique_fly) %>% 
  summarise(genotype = unique(genotype),
            both_wing_index = 100*sum(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ]>(35*pi/180),
                                 na.rm = TRUE)/
              length(Frame[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ]
                     ))
aov1 <- aov(both_wing_index~genotype, data = test_stat_tibble)
TukeyHSD(aov1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = both_wing_index ~ genotype, data = test_stat_tibble)

$genotype
          diff        lwr      upr     p adj
B-A -0.4779248 -2.0714700 1.115620 0.8624926
C-A -0.1148557 -1.4946410 1.264930 0.9963700
D-A  0.8787432 -0.5109769 2.268463 0.3557950
C-B  0.3630691 -1.2026635 1.928802 0.9304240
D-B  1.3566679 -0.2178265 2.931162 0.1170595
D-C  0.9935989 -0.3641396 2.351337 0.2304436

test_stat_tibble <- all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  #filter(genotype != "A") %>% 
  group_by(genotype) %>% 
  group_by(unique_fly) %>% 
  summarise(genotype = unique(genotype),
            both_wing_index = 100*sum(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ]>(15*pi/180),
                                 na.rm = TRUE)/
              length(Frame[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ]
                     ))
aov1 <- aov(both_wing_index~genotype, data = test_stat_tibble)
TukeyHSD(aov1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = both_wing_index ~ genotype, data = test_stat_tibble)

$genotype
          diff       lwr       upr     p adj
B-A  1.3374241 -4.537783  7.212632 0.9338504
C-A  1.9198669 -3.167234  7.006968 0.7588352
D-A 10.3410421  5.217313 15.464771 0.0000040
C-B  0.5824427 -5.190223  6.355108 0.9935904
D-B  9.0036180  3.198649 14.808587 0.0005495
D-C  8.4211753  3.415359 13.426992 0.0001503

all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  group_by(genotype) %>% 
  group_by(unique_fly) %>% 
  summarise(genotype = unique(genotype),
            mean_min_wing = mean(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ])
            ) %>%
  ggplot(aes(x=genotype,y=mean_min_wing)) +
    geom_boxplot() +
    ylim(0,.5)

all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  group_by(genotype) %>% 
  group_by(unique_fly) %>% 
  summarise(genotype = unique(genotype),
            min_wing_gt_15 = sum(mean(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                                                 ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                                        which.max(SmoothedCopulation), 
                                                                                        min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                        ),
                                                                                  min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                                  )
                                             ])>15*pi/180)
            ) %>% 
  summarise(prop_min_wing_gt_15 = 100*sum(min_wing_gt_15>0,na.rm = TRUE)/length(min_wing_gt_15))

all_rawdata %>% 
  filter(dist_to_other__mm > 2) %>% 
  filter(dist_to_other__mm < 10) %>% 
  filter(max_wing_ang__rad > (25*pi/180)) %>%
  #filter(genotype == "D") %>% 
  ggplot(aes(x=genotype, y=wing_l_len__px))+
  geom_boxplot()

Facing angle at start vs end of wing bout

Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'. 
2: Unknown or uninitialised column: 'bin_max_wing'. 
3: Unknown or uninitialised column: 'bin_max_wing'. 
4: Unknown or uninitialised column: 'bin_max_wing'. 
5: Unknown or uninitialised column: 'bin_max_wing'. 
6: Unknown or uninitialised column: 'bin_max_wing'. 
7: Unknown or uninitialised column: 'bin_max_wing'. 
8: Unknown or uninitialised column: 'bin_max_wing'. 

Facing angle at start and end of each bout

temp <- all_rawdata %>% 
  filter(genotype!="CS female") %>% 
  filter(dist_to_other__mm > 2)
  
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
  temp2 <- filter(temp, unique_fly == fly)
  bouts <- rle(temp2$WingGesture)
  
  starts <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
  }
  ends <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    ends[i] = sum(bouts$lengths[1:(i*2)])
  }
  genotype = rep(unique(temp2$genotype),length(which(bouts$values==1)))
  uni_id = rep(unique(temp2$unique_fly),length(which(bouts$values==1)))
  
  temp_tibble <- tibble(unique_fly = uni_id,
                        genotype = genotype,
                        facing_at_start = temp2$facing_angle__rad[starts],
                        facing_at_end = temp2$facing_angle__rad[ends]
                        )
  facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing

Mean facing angle at start and end of bout

temp <- all_rawdata %>% 
  filter(genotype!="CS female") %>% 
  filter(dist_to_other__mm > 2)
  
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
  temp2 <- filter(temp, unique_fly == fly)
  bouts <- rle(temp2$WingGesture)
  
  starts <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
  }
  ends <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    ends[i] = sum(bouts$lengths[1:(i*2)])
  }
  genotype = unique(temp2$genotype)
  uni_id = unique(temp2$unique_fly)
  
  temp_tibble <- tibble(unique_fly = uni_id,
                        genotype = genotype,
                        facing_at_start = mean(temp2$facing_angle__rad[starts]),
                        facing_at_end = mean(temp2$facing_angle__rad[ends])
                        )
  facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing

Median facing angle at start and end of bout (using WingGesture)

temp <- all_rawdata %>% 
  filter(genotype!="CS female") %>% 
  filter(dist_to_other__mm > 2)
  
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
  temp2 <- filter(temp, unique_fly == fly)
  bouts <- rle(temp2$WingGesture)
  
  starts <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
  }
  ends <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    ends[i] = sum(bouts$lengths[1:(i*2)])
  }
  genotype = unique(temp2$genotype)
  uni_id = unique(temp2$unique_fly)
  
  temp_tibble <- tibble(unique_fly = uni_id,
                        genotype = genotype,
                        facing_at_start = median(temp2$facing_angle__rad[starts]),
                        facing_at_end = median(temp2$facing_angle__rad[ends])
                        )
  facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing

Median facing angle at start and end of bout (using >35deg)

temp <- all_rawdata %>% 
  filter(genotype!="CS female") %>% 
  filter(dist_to_other__mm > 2)

Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'. 
2: Unknown or uninitialised column: 'bin_max_wing'. 
3: Unknown or uninitialised column: 'bin_max_wing'. 
4: Unknown or uninitialised column: 'bin_max_wing'. 

  
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
  temp2 <- filter(temp, unique_fly == fly)
  #bouts <- rle(temp2$WingGesture)
  temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad >= (35*pi/180),1,0)
  temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
  bouts <- rle(temp2$bin_max_wing)
  starts <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
  }
  ends <- as.numeric()
  for (i in seq(1,length(which(bouts$values==1)),1)) {
    ends[i] = sum(bouts$lengths[1:(i*2)])
  }
  genotype = unique(temp2$genotype)
  uni_id = unique(temp2$unique_fly)
  
  temp_tibble <- tibble(unique_fly = uni_id,
                        genotype = genotype,
                        facing_at_start = median(temp2$facing_angle__rad[starts]),
                        facing_at_end = median(temp2$facing_angle__rad[ends])
                        )
  facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing

Facing angles with wing extended using WingGesture

Facing angles with wing extended using max_wing>35

test_stat_tibble0 <- all_rawdata %>% 
  filter(genotype!="CS female") %>% 
  group_by(unique_fly) %>% 
  filter(dist_to_other__mm > 2) %>% 
  filter(max_wing_ang__rad > (35*pi/180)) %>%
  summarise(genotype = unique(genotype),
            facing = median(facing_angle__rad)) 
aov0 <- aov(facing~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = facing ~ genotype, data = test_stat_tibble0)

$genotype
            diff        lwr       upr     p adj
B-A -0.064278836 -0.2941883 0.1656306 0.8852971
C-A -0.017473865 -0.2165430 0.1815953 0.9957523
D-A -0.012774271 -0.2132767 0.1877282 0.9983603
C-B  0.046804972 -0.1790918 0.2727017 0.9489797
D-B  0.051504566 -0.1756563 0.2786654 0.9345696
D-C  0.004699594 -0.1911887 0.2005879 0.9999116

test <- all_rawdata %>% 
  filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_1234_2") %>% 
  filter(Id == 36) 

temp2 <- test

Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'. 
2: Unknown or uninitialised column: 'bin_max_wing'. 
3: Unknown or uninitialised column: 'bin_max_wing'. 
4: Unknown or uninitialised column: 'bin_max_wing'. 

#bouts <- rle(temp2$WingGesture)
temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad > (35*pi/180),1,0)
temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
bouts <- rle(temp2$bin_max_wing)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
  starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
  ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = rep(unique(temp2$genotype),length(which(bouts$values==1)))
uni_id = rep(unique(temp2$unique_fly),length(which(bouts$values==1)))
temp_tibble <- tibble(unique_fly = uni_id,
                      genotype = genotype,
                      facing_at_start = temp2$facing_angle__rad[starts],
                      facing_at_end = temp2$facing_angle__rad[ends]
                      )
temp_tibble

starts = for (i in seq(1,length(which(bouts$values==1)),1)) {
                                           starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
                                         }
starts

NULL

test$bin_max_wing <- ifelse(test$max_wing_ang__rad > (35*pi/180),1,0)
test$bin_max_wing <- test$bin_max_wing %>% replace_na(0)
#test$bin_max_wing
bouts <- rle(test$bin_max_wing)
bouts

Run Length Encoding
  lengths: int [1:393] 675 3 2 14 2 22 212 15 15 2 ...
  values : num [1:393] 0 1 0 1 0 1 0 1 0 1 ...

bouts <- rle(test$WingGesture)
bouts

Run Length Encoding
  lengths: int [1:1361] 666 5 3 45 198 1 2 5 4 17 ...
  values : num [1:1361] 0 1 0 1 0 1 0 1 0 1 ...

sum(test$bin_max_wing[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)]==test$WingGesture[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)])

[1] 1230

sum(test$bin_max_wing[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)]!=test$WingGesture[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)])

[1] 35

gg_color_hue(4)

[1] "#F8766D" "#7CAE00" "#00BFC4" "#C77CFF"

max(test$Frame)

[1] 5974

temp2 <- all_rawdata %>% 
  filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_123_3") %>% 
  filter(Id == 1) %>% 
  filter(dist_to_other__mm > 2)
#bouts <- rle(temp2$WingGesture)
temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad >= (35*pi/180),1,0)
temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
bouts <- rle(test$bin_max_wing)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
  starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
  ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = unique(temp2$genotype)
uni_id = unique(temp2$unique_fly)
temp_tibble <- tibble(unique_fly = uni_id,
                      genotype = genotype,
                      facing_at_start = median(temp2$facing_angle__rad[starts]),
                      facing_at_end = median(temp2$facing_angle__rad[ends])
                      )
temp_tibble

bouts <- rle(test$bin_max_wing)
bouts

Run Length Encoding
  lengths: int [1:198] 631 3 57 5 363 34 357 1 23 1 ...
  values : num [1:198] 0 1 0 1 0 1 0 1 0 1 ...

Trying to make functions for courtship window and calculate indices… they aren’t working right, and seem to be ignoring the grouping of the tibble and returning the same value for every individual…

If i have to for loop over the groups of a tibble, the below line might be usefull. dim(unique(all_rawdata[all_rawdata %>% group_by(genotype) %>% group_vars()]))[1]

courtship_window <- function(input,wind=600,...){
  temp_tibble <- summarise(input, start_of_courtship = which.max(SmoothedCourtship),
                                  end_of_courtship = ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                             ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                    which.max(SmoothedCopulation), 
                                                                    min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                    ),
                                                              min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                              )
                        )
  return(temp_tibble$start_of_courtship:temp_tibble$end_of_courtship)
  #return(c(temp_tibble$start_of_courtship,temp_tibble$end_of_courtship))
  #return(temp_tibble)
}

courtship_window <- function(input,wind=600,...){
  temp_tibble <- summarise(input, start_of_courtship = which.max(SmoothedCourtship),
                                  end_of_courtship = ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                             ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                    which.max(SmoothedCopulation), 
                                                                    min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                    ),
                                                              min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                              )
                        )
  #return(temp_tibble$start_of_courtship:temp_tibble$end_of_courtship)
  return(c(temp_tibble$start_of_courtship,temp_tibble$end_of_courtship))
  #return(temp_tibble)
}

courtship_window(test, wind = 600)

temp <- all_rawdata %>% 
  filter(genotype!="CS female")

uniq_fly <- unique(temp$unique_fly)
courting_frames <- tibble()
for (fly in uniq_fly) {
  temp2 <- temp %>% filter(unique_fly == fly) 
  temp2 <- slice(temp2, which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
                                                             ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)), 
                                                                    which.max(SmoothedCopulation), 
                                                                    min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                                    ),
                                                              min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
                                                              ))
  courting_frames <- bind_rows(courting_frames,temp2)
}
courting_frames

calculate_indices <- function(input = .,feature,jaaba=TRUE,thresh=NULL,wind=600){
  if (feature %in% names(input)) {
    var_ind = paste0(feature,"_index")
    if (jaaba) {
      temp_tibble <- summarise(input, genotype = unique(genotype), 
                                      !!var_ind := 100*sum(input[[feature]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
                                                       ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)), 
                                                              which.max(input[["SmoothedCopulation"]]), 
                                                              min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                              ),
                                                        min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                        )])/
                                                    length(input[["Frame"]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
                                                       ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)), 
                                                              which.max(input[["SmoothedCopulation"]]), 
                                                              min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                              ),
                                                        min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                        )])
                )
    } else {
      temp_tibble <- summarise(input, genotype = unique(genotype), 
                                      !!var_ind := 100*sum(input[[feature]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
                                                       ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)), 
                                                              which.max(input[["SmoothedCopulation"]]), 
                                                              min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                              ),
                                                        min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                        )]>thresh,na.rm = TRUE)/
                                                    length(input[["Frame"]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
                                                       ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)), 
                                                              which.max(input[["SmoothedCopulation"]]), 
                                                              min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                              ),
                                                        min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
                                                        )])
                )
    }
    return(temp_tibble)
  } else {
    message(paste0(feature, " does not exist in table"))
  }
} 

calculate_indices <- function(input,feature,jaaba=TRUE,thresh=NULL){
  if (feature %in% names(input)) {
    var_ind = paste0(feature,"_index")
    if (jaaba) {
      temp_tibble <- summarise(input, genotype = unique(genotype), 
                                      !!var_ind := 100*sum(input[[feature]][courtship_window(input)])/
                                                    length(input[["Frame"]][courtship_window(input)])
                )
    } else {
      temp_tibble <- summarise(input, genotype = unique(genotype), 
                                      !!var_ind := 100*sum(input[[feature]][courtship_window(input)]>thresh,na.rm = TRUE)/
                                                    length(input[["Frame"]][courtship_window(input)])
                )
    }
    return(temp_tibble)
  } else {
    message(paste0(feature, " does not exist in table"))
  }
} 

all_rawdata %>% 
  filter(genotype != "CS female") %>% 
  filter(genotype == "D") %>% 
  group_by(genotype) %>% 
  group_by(unique_fly) %>% 
  calculate_indices(feature = "max_wing_ang__rad",jaaba = FALSE,thresh = (25*pi/180))

test <- all_rawdata %>% 
  filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_1234_2") %>% 
  filter(Id == 36) 

calculate_indices(input = test,feature = "WingGesture",jaaba = TRUE)

calculate_indices(input = test,feature = "max_wing_ang__rad",jaaba = FALSE,thresh = (25*pi/180))

calculate_indices(input = test,feature = "max_wing_angle__rad",jaaba = FALSE,thresh = (25*pi/180))